import numpy as np
from scipy.integrate import odeint
import pandas as pd
from sko.SA import SA
import matplotlib.pyplot as plt

# 1. 水的粘度查找表 (单位: mPa*s)
water_viscosity = {
    30: 0.8007, 31: 0.7840, 32: 0.7679, 33: 0.7523, 34: 0.7371,
    35: 0.7225, 36: 0.7085, 37: 0.6947, 38: 0.6814, 39: 0.6685,
    40: 0.6560, 41: 0.6439, 42: 0.6321, 43: 0.6207, 44: 0.6097,
    45: 0.5988, 46: 0.5883, 47: 0.5782, 48: 0.5683, 49: 0.5588,
    50: 0.5494
}


# 2. 单位转换函数
def convert_viscosity(vis):
    """mPa*s → Pa*s"""
    return vis * 0.001


# 3. 参数配置
k_params = {
    'k1': 2.12911238, 'kH': 0.74896823, 'k3': 0.1803147,
    'k4': 0.56862499, 'k5': 4.68784022, 'k6': 0.56668138,
    'S_C': 0.3878337, 'S_S': 0.64657911, 'alpha': 0.31177578
}
viscosity_rate = 0.96  # DMF/水粘度比
ro_D = 948  # kg/m³
ro_S = 1261  # kg/m³
ro_C = 1300  # kg/m³
m_D0 = 24e-3  # kg
m_S = 6e-3  # kg
r = 0.3413e-9  # m (环丁砜分子半径)


# 4. 目标函数
def calculate_pore_area(x):
    T, H, solid_content = x
    print(f"温度:{T}, 湿度:{H}, 固含量:{solid_content}\r\n")
    if T > 50 or T < 30 or H > 90 or T < 50 or solid_content > 0.1 or T < 0.06:
        return 0
    # 计算醋酸纤维素初始质量
    m_C = (m_D0 + m_S) * solid_content

    # 计算总体积 (m³)
    V = m_D0 / ro_D + m_S / ro_S + m_C / ro_C

    # 计算饱和蒸气压 (Pa)
    T_kelvin = T + 273.15
    Psat = 10 ** (6.09451 - 2725.96 / (T_kelvin + 28.209)) * 0.133322

    # 计算蒸发系数 A0
    A0 = -k_params['k1'] * Psat * (1 - k_params['kH'] * H / 100) / ro_D / V

    # 终止时间（当DMF质量降至alpha倍时）
    end_time = np.log(k_params['alpha']) / A0

    # 时间网格 (自适应步长)
    if end_time > 0.01:
        t = np.linspace(0, end_time, min(1000, int(end_time * 100)))
    else:
        return 0  # 避免无效计算

    # DMF质量函数
    m_D_func = lambda t: m_D0 * np.exp(A0 * t)

    # 析出质量计算
    m_S_out = lambda t: max(0, m_S - m_D_func(t) * k_params['S_S'])
    m_C_out = lambda t: max(0, m_C - m_D_func(t) * k_params['S_C'])
    phi_C_out = lambda t: m_C_out(t) / (ro_C * V)

    # 在eta计算中使用：
    T_int = int(round(T))
    eta0 = convert_viscosity(water_viscosity[T_int]) * viscosity_rate

    # 溶液粘度
    eta = lambda t: eta0 * (1 + 2.5 * phi_C_out(t))

    # 扩散系数
    r = 0.3413e-9  # 环丁砜分子半径 (m)
    D = lambda t: 1.38e-23 * T_kelvin / (6 * np.pi * r * eta(t))

    # 液滴速度
    v = lambda t: k_params['k3'] * np.sqrt(D(t))

    # 分子数密度
    n_density = lambda t: m_S_out(t) / (ro_S * (4 / 3) * np.pi * r ** 3) / V

    # 碰撞频率
    Z = lambda t: np.sqrt(2) * n_density(t) * np.pi * r ** 2 * v(t)

    # 微分方程定义
    def dmdt(m_s, t):
        return (k_params['k4'] * Z(t) +
                k_params['k5'] * v(t) * (m_S_out(t) - m_s) / V * m_s)

    # 求解微分方程
    sol = odeint(dmdt, 0, t)
    pore_area = k_params['k6'] * sol[-1, 0]
    print(f"温度:{T}, 湿度:{H}, 固含量:{solid_content}, 最终结果:{pore_area}\r\n")
    # 确保所有数学运算使用浮点数
    return -float(pore_area)  # 负值用于最小化


# 初始化模拟退火算法
sa = SA(
    func=calculate_pore_area,      # 目标函数
    x0=[40, 70, 0.08],        # 初始解（随机起点）
    T_max=1,             # 初始温度（高温）
    T_min=1e-9,          # 终止温度（低温）
    L=300,               # 每个温度下的迭代次数
    max_stay_counter=150 # 停止条件：若连续多次无改进则终止
)

# 3. 运行算法
best_x, best_y = sa.run()

# 4. 输出结果
print('best_x:', best_x)
print('best_y:', best_y)

# 5. 可视化优化过程（历史最优值曲线）
plt.plot(pd.DataFrame(sa.best_y_history).cummin(axis=0))
plt.xlabel('iter')
plt.ylabel('target')
plt.title('SA')
plt.show()
